Deaths by horse kicks
Ladislaus Josephovich Bortkiewicz was a Russian economist and statistician. He noted that the Poisson distribution can be very useful in applied statistics when describing low-frequency events in a large population. In a well-known example, he showed that the number of deaths by horse kicks in the Prussian army followed the Poisson distribution.
Consider the following two sets of observations, taken over a fixed large time interval in two different corps:
| \(y\) dead soldiers | 0 | 1 | 2 | 3 | 4 | ≥ 5 |
|---|---|---|---|---|---|---|
| \(n_1\) observations | 109 | 65 | 22 | 3 | 1 | 0 |
| \(n_2\) observations | 144 | 91 | 32 | 11 | 2 | 0 |
The two datasets are independent, so we can safely join them. Another option could be to analyse one before the other, taking the posterior from the first as prior for the second: but since we’ll be using Gamma priors, the update rule of the prior parameters,
\[ \alpha_{\mathrm{post}} = \alpha_{\mathrm{prior}} + \sum_i y_i, \quad \lambda_{\mathrm{post}} = \alpha_{\mathrm{prior}} + n \]
would give the same result as for the case with the two datasets lumped together.
We’ll assume a uniform prior, and compute the posterior distribution for \(\lambda\), the death rate over the measurement time. We’ll do the same also with a Jeffrey’s prior.
y <- c(
rep(0, 109), rep(1, 65), rep(2, 22), rep(3, 3), rep(4, 1),
rep(0, 144), rep(1, 91), rep(2, 32), rep(3, 11), rep(4, 2)
)
# define a common function to return the Bayes stuff
# s = sum of observations
# n = number of observations
# pr_* = prior parameters
baypois <- function(s, n, pr_alpha, pr_lambda) {
alpha <- pr_alpha + s
lambda <- pr_lambda + n
post <- function(x) dgamma(x, alpha, lambda)
median <- qgamma(0.5, alpha, lambda)
mean <- integrate(\(x) x * post(x), 0, Inf)$value
std <- sqrt(integrate(\(x) (x - mean)^2 * post(x), 0, Inf)$value)
cred <- sapply(c(0.025, 0.975), \(c) qgamma(c, alpha, lambda))
list(post = post, med = median, mean = mean, std = std, cred = cred)
}
unif <- baypois(sum(y), length(y), 1, 0)
jeff <- baypois(sum(y), length(y), 0.5, 0)| Prior | Median | Mean | St. dev. | 95% cred. int. |
|---|---|---|---|---|
| Uniform | 0.664 | 0.665 | 0.0372 | 0.594, 0.739 |
| Jeffrey’s | 0.663 | 0.664 | 0.0372 | 0.593, 0.738 |
The parameters of the two posteriors are almost identical, but we could have expected it given the high number of observations (480).
Since the distribution are basically overlapping, we’ll zoom on the region
my_pal <- wesanderson::wes_palette("Zissou1", 5)[c(1, 3, 5)]
posts <- tibble(
lam = seq(0.5, 0.9, by = 0.0005),
unif = unif$post(lam), jeff = jeff$post(lam)
) |>
pivot_longer(-lam, names_to = "dist", values_to = "prob") |>
mutate(dist = fct_relevel(dist, "unif", "jeff"))
modes <- posts |>
group_by(dist) |>
summarize(mode = lam[which.max(prob)], pmax = max(prob))
creds <- posts |>
group_by(dist) |>
filter(
lam > lam[which.max(cumsum(prob) * 0.0005 > 0.025)] &
lam < lam[which.max(cumsum(prob) * 0.0005 > 0.975)]
)
posts |>
ggplot(aes(x = lam, y = prob)) +
geom_line(linewidth = 0.8, colour = my_pal[1]) +
geom_area(
data = creds,
fill = my_pal[1], alpha = 0.5,
position = "identity", show.legend = FALSE
) +
geom_segment(
aes(x = mode, y = 0, xend = mode, yend = pmax),
data = modes, colour = my_pal[3], linewidth = 0.8
) +
scale_colour_manual(
values = my_pal[c(1, 3)], labels = c("Uniform", "Jeffrey’s")
) +
geom_label(
aes(
x = mode, y = 0,
label = paste("Mode =", format(mode, digits = 3))
),
data = modes, hjust = -0.1, vjust = -0.3, colour = my_pal[3]
) +
labs(
x = "Poisson parameter λ",
y = "Posterior probability P(λ | y)",
) +
facet_wrap(
vars(dist), nrow = 2,
labeller = as_labeller(c(
unif = "With uniform prior",
jeff = "With Jeffrey’s prior"
))
)
Deaths by horse kicks, but with JAGS
Now we’ll repeat the previous analysis, but this time using a Markov
Chain Monte Carlo sampling. We’ll rely on the JAGS library
– and its R interface rjags – to build the chain and
generate samples with the Gibbs method.
library(rjags)
model <- "model {
# data likelihood
for (i in 1:length(x)) {
x[i] ~ dpois(lambda)
}
# uniform prior for lambda
lambda ~ dexp(0.0001)
# predicted data, given lambda
y ~ dpois(lambda)
}"
data <- list(x = c(
rep(0, 109), rep(1, 65), rep(2, 22), rep(3, 3), rep(4, 1),
rep(0, 144), rep(1, 91), rep(2, 32), rep(3, 11), rep(4, 2)
))
# create the model
jm <- jags.model(textConnection(model), data)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 480
## Unobserved stochastic nodes: 2
## Total graph size: 483
##
## Initializing model# update the Markov chain (burn-in)
update(jm, 1000)
chain <- coda.samples(jm, c("lambda", "y"), n.iter = 1e5)
plot(chain, col = "turquoise4", right = FALSE)## Warning in plot.window(...): "right" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "right" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "right" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "right" is not a
## graphical parameter## Warning in box(...): "right" is not a graphical parameter## Warning in title(...): "right" is not a graphical parameter## Warning in plot.window(...): "right" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "right" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "right" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "right" is not a
## graphical parameter## Warning in box(...): "right" is not a graphical parameter## Warning in title(...): "right" is not a graphical parameter